package com.ruoyi.huanjing.tool;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.interpolation.BicubicSplineInterpolator;
import org.apache.commons.math3.analysis.interpolation.BicubicSplineInterpolatingFunction;
import org.apache.commons.math3.util.MathUtils;

public class DataSmoothingUtils {

    /**
     * 高斯滤波平滑二维数据
     * @param data 原始二维数据
     * @param sigmaTime 时间维度高斯核标准差
     * @param sigmaHeight 高度维度高斯核标准差
     * @return 平滑后的数据
     */
    public static double[][] gaussianFilter(double[][] data, double sigmaTime, double sigmaHeight) {
        int rows = data.length;
        if (rows == 0) return data;
        int cols = data[0].length;
        double[][] result = new double[rows][cols];

        // 生成高斯核
        int kernelSizeTime = Math.max(3, (int) (6 * sigmaTime + 1));
        int kernelSizeHeight = Math.max(3, (int) (6 * sigmaHeight + 1));
        double[][] kernelTime = generateGaussianKernel(kernelSizeTime, sigmaTime);
        double[][] kernelHeight = generateGaussianKernel(kernelSizeHeight, sigmaHeight);

        // 先对时间维度卷积
        double[][] temp = new double[rows][cols];
        for (int i = 0; i < rows; i++) {
            temp[i] = convolve1D(data[i], kernelTime);
        }

        // 再对高度维度卷积
        for (int j = 0; j < cols; j++) {
            double[] colData = new double[rows];
            for (int i = 0; i < rows; i++) {
                colData[i] = temp[i][j];
            }
            double[] smoothedCol = convolve1D(colData, kernelHeight);
            for (int i = 0; i < rows; i++) {
                result[i][j] = smoothedCol[i];
            }
        }

        return result;
    }

    private static double[][] generateGaussianKernel(int size, double sigma) {
        double[][] kernel = new double[size][size];
        double sum = 0.0;
        int center = size / 2;
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < size; j++) {
                double x = i - center;
                double y = j - center;
                kernel[i][j] = Math.exp(-(x * x + y * y) / (2 * sigma * sigma)) 
                    / (2 * Math.PI * sigma * sigma);
                sum += kernel[i][j];
            }
        }
        // 归一化
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < size; j++) {
                kernel[i][j] /= sum;
            }
        }
        return kernel;
    }

    private static double[] convolve1D(double[] data, double[][] kernel) {
        int n = data.length;
        int m = kernel.length;
        int half = m / 2;
        double[] result = new double[n];
        for (int i = 0; i < n; i++) {
            double sum = 0.0;
            for (int j = 0; j < m; j++) {
                int k = i + j - half;
                if (k >= 0 && k < n) {
                    sum += data[k] * kernel[j][j];
                }
            }
            result[i] = sum;
        }
        return result;
    }

    /**
     * 双三次插值平滑数据（增加分辨率）
     * @param data 原始二维数据
     * @param originalTime 原始时间轴数据
     * @param originalHeight 原始高度轴数据
     * @param newTimeLength 新的时间轴长度
     * @param newHeightLength 新的高度轴长度
     * @return 包含平滑后数据、新时间轴、新高度轴的对象数组
     */
    public static Object[] bicubicInterpolation(double[][] data, double[] originalTime, double[] originalHeight,
                                                int newTimeLength, int newHeightLength) {
        // 验证输入数据
        if (data == null || data.length == 0 || originalTime == null || originalHeight == null) {
            throw new IllegalArgumentException("输入数据不能为空");
        }
        if (data.length != originalTime.length || data[0].length != originalHeight.length) {
            throw new IllegalArgumentException("数据维度与轴数据不匹配");
        }

        BicubicSplineInterpolator interpolator = new BicubicSplineInterpolator();
        // 关键修正：使用正确的返回类型
        BicubicSplineInterpolatingFunction interpolatingFunction = 
            interpolator.interpolate(originalTime, originalHeight, data);

        // 生成新的时间轴
        double[] newTime = new double[newTimeLength];
        double timeMin = originalTime[0];
        double timeMax = originalTime[originalTime.length - 1];
        for (int i = 0; i < newTimeLength; i++) {
            newTime[i] = timeMin + (timeMax - timeMin) * i / (newTimeLength - 1);
        }

        // 生成新的高度轴
        double[] newHeight = new double[newHeightLength];
        double heightMin = originalHeight[0];
        double heightMax = originalHeight[originalHeight.length - 1];
        for (int i = 0; i < newHeightLength; i++) {
            newHeight[i] = heightMin + (heightMax - heightMin) * i / (newHeightLength - 1);
        }

        // 计算插值后的平滑数据
        double[][] smoothedData = new double[newTimeLength][newHeightLength];
        for (int i = 0; i < newTimeLength; i++) {
            for (int j = 0; j < newHeightLength; j++) {
                // 关键修正：使用正确的方法获取插值结果
                smoothedData[i][j] = interpolatingFunction.value(newTime[i], newHeight[j]);
            }
        }

        return new Object[]{smoothedData, newTime, newHeight};
    }
}
    